Models

In this document, we will fit all models regarding this project. The detailed preregistered analysis plan can be found in our OSF’s project.

Data

Let’s load libraries and data.

# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")

# Install/Load packages

pacman::p_load("rio",
               "here",
               "metafor",
               "dplyr",
               "bayesmeta",
               "brms")

remotes::install_github("stan-dev/cmdstanr")

# Ensure that packages are installed according to versions in renv's lockfile
renv::restore()

# Import data
d = rio::import(here::here("data", "data.xlsx"))

Let’s calculate the crude odds ratio based on the extracted data.

# Calculate crude log odds ratios

d_logOR = 
  metafor::escalc(
  measure = "OR", # log odds ratio,
  
  # Tocilizumab
  ai = trt_events,
  n1i = trt_total,
  
  # Control
  ci = control_events,
  n2i = control_total,
  
  data = d
) %>%
  as_tibble() %>% 
  # Calculate standard error
  dplyr::mutate(sei = sqrt(vi),
         # Set order to use "tocilizumab" as the Intercept in brms
         treatment = factor(treatment,
                         levels = c("tocilizumab",
                                    "sarilumab",
                                    "baricitinib")))

# saveRDS(d_logOR,
#         here("output", "data", "effect_sizes.Rds"))

Meta-analyses

First, we will fit a single meta-analysis model per drug (tocilizumab, sarilumab, or baricitinib vs. control).

All Bayesian meta-analysis random-effect models are the same, and is described as:

\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood} \\ \theta_i & \sim Normal(\mu, \tau^2)\\ \\ \mu & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors} \\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2) \\ \end{align*} \]

where \(y_i\) is the observed mean log odds ratio in study \(i\) for either tocilizumab, sarilumab, or baricitinib versus control treatment. We assume these effect sizes are normally distributed around the study-specific mean \(\theta_i\) along with a known sampling variance, represented by the observed \(\sigma_i^2\). We also assume \(\theta_i\) is drawn from a normal distribution where \(\mu\) is the average effect and \(\tau^2\) is the between-study heterogeneity.

## Formula
mf = 
  # https://bookdown.org/content/4857/horoscopes-insights.html#consider-using-the-0-intercept-syntax
  formula(yi | se(sei) ~ 0 + Intercept + (1 | study))

## Priors

# Tau
informative = bayesmeta::TurnerEtAlPrior("all-cause mortality",
                                         "pharmacological",
                                         "placebo / control")


logmean = informative$parameters["tau", "mu"]
logsd = informative$parameters["tau", "sigma"]

# logmean
# -1.975

# logsd
# 0.67

priors_ma = 
  brms::prior(normal(0, 0.75), class = "b", coef = "Intercept") +
  brms::prior(lognormal(-1.975, 0.67), class = "sd") 

## Models

# Tocilizumab

ma_toci = 
  brms::brm(
  data = d_logOR |> dplyr::filter(treatment == "tocilizumab"),
  family = gaussian,
  
  formula = mf,
  prior = priors_ma,
  sample_prior = TRUE,
  
  control = list(adapt_delta = .95),
  backend = "cmdstanr", # faster
  cores = parallel::detectCores(),
  chains = 4,
  warmup = 2000, 
  iter = 4000, 
  seed = 123,
  
  file = here::here("output", "fits", "ma_toci.Rds"),
  file_refit = "on_change"
)

# Sarilumab

ma_sari = 
  brms::brm(
  data = d_logOR |> dplyr::filter(treatment == "sarilumab"),
  family = gaussian,
  
  formula = mf,
  prior = priors_ma,
  sample_prior = TRUE,
  
  control = list(adapt_delta = .95),
  backend = "cmdstanr", # faster
  cores = parallel::detectCores(),
  chains = 4,
  warmup = 2000, 
  iter = 4000, 
  seed = 123,
  
  file = here::here("output", "fits", "ma_sari.Rds"),
  file_refit = "on_change"
)

# Baricitinib

ma_bari = 
  brms::brm(
  data = d_logOR |> dplyr::filter(treatment == "baricitinib", study != "RECOVERY Bari"),
  family = gaussian,
  
  formula = mf,
  prior = priors_ma,
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  cores = parallel::detectCores(),
  chains = 4,
  warmup = 2000, 
  iter = 4000, 
  seed = 123,
  
  file = here::here("output", "fits", "ma_bari.Rds"),
  file_refit = "on_change"
)

Sensitivity analysis with all patients in RECOVERY Bari

ma_bari_sens = 
  brms::brm(
  data = d_logOR |> dplyr::filter(treatment == "baricitinib", study != "RECOVERY Bari (No Toci)"),
  family = gaussian,
  
  formula = mf,
  prior = priors_ma,
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  cores = parallel::detectCores(),
  chains = 4,
  warmup = 2000, 
  iter = 4000, 
  seed = 123,
  
  file = here::here("output", "fits", "ma_bari_sens.Rds"),
  file_refit = "on_change"
)

Meta-regressions

Now, we will fit Bayesian meta-regression models (one for tocilizumab vs. sarilumab data, another for tocilizumab vs. baricitinib).

The model can be described as:

\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood}\\ \theta_i & \sim Normal(\mu, \tau^2)\\ \mu &= \beta_0 + \beta_1 x_i\\ \\ \end{align*} \]

where \(y_i\) is the observed mean log odds ratio in study \(i\) for either tocilizumab, sarilumab, or baricitinib versus control treatment. We assume these effect sizes are normally distributed around the study-specific mean \(\theta_i\) along with a known sampling variance, represented by the observed \(\sigma_i^2\). We also assume \(\theta_i\) is drawn from a normal distribution where \(\mu\) is the average effect and \(\tau^2\) is the between-study heterogeneity. We will estimate \(\mu\) as a function of the estimated population effect size \(\beta_0\) and the moderator \(x\) multiplied by the coefficient \(\beta_1\). \(x\) is dummy-coded, where indicates that the observed \(y_i\) is an effect size estimate of compared to control, while indicates that \(y_i\) is an effect size estimate of the compared to control.

Tocilizumab vs. Sarilumab

Here are the priors for this comparison:

\[ \begin{align*} \beta_0 & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors}\\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2)\\ \beta_{1[Skeptical]} & \sim \operatorname{Normal}(0, 0.354^2)\\ \beta_{1[Sarilumab]} & \sim \operatorname{Normal}(-0.049, 0.193^2)\\ \beta_{1[Tocilizumab]} & \sim \operatorname{Normal}(0.049, 0.193^2)\\ \beta_{1[Vague]} & \sim \operatorname{Normal}(0, 4^2) \end{align*} \]

## Means

# Figure S3 in https://doi.org/10.1101/2021.06.18.21259133;
# version posted June 22, 2021
reciprocal_remap_cap = 1/1.05
mean_log_sari = log(reciprocal_remap_cap)

mean_skeptical = log(1)

## SDs

# Assuming a mean = log(1), what is the SD that yields a distribution
# with 0.025 probability density below log(0.5)?
sari_sd_skeptical = (mean_skeptical - log(0.5))/qnorm(1-0.025)

# Assuming a mean = log(0.952), what is the SD that yields a distribution
# with 0.4 probability density above log(1)?
sari_sd_optimistic_sari = (log(1) - mean_log_sari)/qnorm(1-0.4)

# Assuming a mean = log(1/0.952) = -log(0.952), what is the SD that yields a
# distribution with 0.4 probability density below log(1)?
sari_sd_optimistic_toci = (-mean_log_sari - log(1))/qnorm(1-0.4)

sari_sd_vague = 4

# General model characteristics 

## Formula
mf = 
  formula(yi | se(sei) ~ 0 + Intercept + treatment + (1 | study))

## Priors for beta_0 and tau

# tau:
# logmean
# -1.975

# logsd
# 0.67

priors = 
  brms::prior(normal(0, 0.75), class = "b", coef = "Intercept") +
  brms::prior(lognormal(-1.975, 0.67), class = "sd", group = "study")

Fit models

## Skeptical model

# sari_sd_skeptical
# 0.353653

m_sari_skeptical =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentsarilumab"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_sari_skeptical.Rds"),
  file_refit = "on_change"
)

## "Optimistic for Sarilumab" model

# mean_log_sari
# -0.04879016

# sari_sd_optimistic_sari
# 0.1925823

m_sari_optimistic_sari =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(-0.04879016, 0.1925823), class = "b", coef = "treatmentsarilumab"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_sari_optimistic_sari.Rds"),
  file_refit = "on_change"
)

## "Optimistic for Tocilizumab" model

# -mean_log_sari
# 0.04879016

# sari_sd_optimistic_toci
# 0.1925823

m_sari_optimistic_toci =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0.04879016, 0.1925823), class = "b", coef = "treatmentsarilumab"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_sari_optimistic_toci.Rds"),
  file_refit = "on_change"
)

## Vague model

m_sari_vague =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 4), class = "b", coef = "treatmentsarilumab"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_sari_vague.Rds"),
  file_refit = "on_change"
)

Tocilizumab vs. Baricitinib

Here are the priors for this comparison:

\[ \begin{align*} \beta_0 & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors}\\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2)\\ \beta_{1[Skeptical]} & \sim \operatorname{Normal}(0, 0.354^2)\\ \beta_{1[Baricitinib]} & \sim \operatorname{Normal}(-0.105, 0.416^2)\\ \beta_{1[Tocilizumab]} & \sim \operatorname{Normal}(0.105, 0.416^2)\\ \beta_{1[Vague]} & \sim \operatorname{Normal}(0, 4^2) \end{align*} \]

## Means
mean_log_bari = log(0.9)

mean_skeptical = log(1)

## SDs

# Assuming a mean = log(1), what is the SD that yields a distribution
# with 0.025 probability density below log(0.5)?
bari_sd_skeptical = (mean_skeptical - log(0.5))/qnorm(1-0.025)

# Assuming a mean = log(0.9), what is the SD that yields a distribution
# with 0.4 probability density above log(1)?
bari_sd_optimistic_bari = (log(1) - mean_log_bari)/qnorm(1-0.4)

# Assuming a mean = log(1/0.9) = -log(0.9), what is the SD that yields a
# distribution with 0.4 probability density below log(1)?
bari_sd_optimistic_toci = (-mean_log_bari - log(1))/qnorm(1-0.4)

bari_sd_vague = 4

## Skeptical model

# bari_sd_skeptical
# 0.353653

m_bari_skeptical =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_skeptical.Rds"),
  file_refit = "on_change"
)


## "Optimistic for Baricitinib" model

# mean_log_bari
# -0.1053605

# bari_sd_optimistic_bari
# 0.4158742

m_bari_optimistic_bari =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(-0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_optimistic_bari.Rds"),
  file_refit = "on_change"
)

## "Optimistic for Tocilizumab" model

# -mean_log_bari
# 0.1053605

# bari_sd_optimistic_toci
# 0.4158742

m_bari_optimistic_toci =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_optimistic_toci.Rds"),
  file_refit = "on_change"
)

## Vague model

m_bari_vague =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 4), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_vague.Rds"),
  file_refit = "on_change"
)

Sensitivity analysis with all patients in RECOVERY Bari

## Skeptical model

# bari_sd_skeptical
# 0.353653

m_bari_skeptical_sens =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_skeptical_sens.Rds"),
  file_refit = "on_change"
)


## "Optimistic for Baricitinib" model

# mean_log_bari
# -0.1053605

# bari_sd_optimistic_bari
# 0.4158742

m_bari_optimistic_bari_sens =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(-0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_optimistic_bari_sens.Rds"),
  file_refit = "on_change"
)

## "Optimistic for Tocilizumab" model

# -mean_log_bari
# 0.1053605

# bari_sd_optimistic_toci
# 0.4158742

m_bari_optimistic_toci_sens =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_optimistic_toci_sens.Rds"),
  file_refit = "on_change"
)

## Vague model

m_bari_vague_sens =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 4), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_vague_sens.Rds"),
  file_refit = "on_change"
)

Model diagnostics

# Load function
source(here::here("functions", "diag_plot.R"))

Meta-analyses

Tocilizumab model

diag_plot(model = ma_toci,
          pars_list = c("b_Intercept", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(ma_toci,
               type = "dens_overlay",
               ndraws = 50)

brms::pp_check(ma_toci,
               type = "ecdf_overlay",
               ndraws = 50)

Sarilumab model

diag_plot(model = ma_sari,
          pars_list = c("b_Intercept", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(ma_sari,
               type = "dens_overlay",
               ndraws = 50)

brms::pp_check(ma_sari,
               type = "ecdf_overlay",
               ndraws = 50)

Baricitinib models

diag_plot(model = ma_bari,
          pars_list = c("b_Intercept", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(ma_bari,
               type = "dens_overlay",
               ndraws = 50) + coord_cartesian(y = c(0, 5))

brms::pp_check(ma_bari,
               type = "ecdf_overlay",
               ndraws = 50)

Baricitinib models

diag_plot(model = ma_bari_sens,
          pars_list = c("b_Intercept", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(ma_bari_sens,
               type = "dens_overlay",
               ndraws = 50) + coord_cartesian(y = c(0, 5))

brms::pp_check(ma_bari_sens,
               type = "ecdf_overlay",
               ndraws = 50)

Meta-regressions

Tocilizumab vs. Sarilumab

“Skeptical” model

diag_plot(model = m_sari_skeptical,
          pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_sari_skeptical,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50)

brms::pp_check(m_sari_skeptical,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Optimistic for Sarilumab” model

diag_plot(model = m_sari_optimistic_sari,
          pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_sari_optimistic_sari,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50)

brms::pp_check(m_sari_optimistic_sari,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Optimistic for Tocilizumab” model

diag_plot(model = m_sari_optimistic_toci,
          pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_sari_optimistic_toci,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50)

brms::pp_check(m_sari_optimistic_toci,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Vague” model

diag_plot(model = m_sari_vague,
          pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_sari_vague,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50)

brms::pp_check(m_sari_vague,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

Tocilizumab vs. Baricitinib

“Skeptical” models

diag_plot(model = m_bari_skeptical,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_skeptical,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_skeptical,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

diag_plot(model = m_bari_skeptical_sens,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_skeptical_sens,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_skeptical_sens,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Optimistic for Baricitinib” models

diag_plot(model = m_bari_optimistic_bari,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_optimistic_bari,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_optimistic_bari,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

diag_plot(model = m_bari_optimistic_bari_sens,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_optimistic_bari_sens,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_optimistic_bari_sens,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Optimistic for Tocilizumab” models

diag_plot(model = m_bari_optimistic_toci,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_optimistic_toci,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_optimistic_toci,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Optimistic for Tocilizumab” model

diag_plot(model = m_bari_optimistic_toci_sens,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_optimistic_toci_sens,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_optimistic_toci_sens,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Vague” models

diag_plot(model = m_bari_vague,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_vague,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_vague,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

diag_plot(model = m_bari_vague_sens,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_vague_sens,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_vague_sens,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)